home *** CD-ROM | disk | FTP | other *** search
/ AGA Toolkit '97 / The AGA Toolkit '97.iso / miscellaneous / science / maths / calc / source / qfunc.c < prev    next >
Encoding:
C/C++ Source or Header  |  1996-09-07  |  24.0 KB  |  1,219 lines
  1. /*
  2.  * Copyright (c) 1994 David I. Bell
  3.  * Permission is granted to use, distribute, or modify this source,
  4.  * provided that this copyright notice remains intact.
  5.  *
  6.  * Extended precision rational arithmetic non-primitive functions
  7.  */
  8.  
  9. #include "qmath.h"
  10.  
  11.  
  12. NUMBER *_epsilon_;    /* default precision for real functions */
  13. long _epsilonprec_;    /* bits of precision for epsilon */
  14.  
  15.  
  16. /*
  17.  * Set the default precision for real calculations.
  18.  * The precision must be between zero and one.
  19.  */
  20. void
  21. setepsilon(q)
  22.     NUMBER *q;        /* number to be set as the new epsilon */
  23. {
  24.     NUMBER *old;
  25.  
  26.     if (qisneg(q) || qiszero(q) || (qreli(q, 1L) >= 0))
  27.         math_error("Epsilon value must be between zero and one");
  28.     old = _epsilon_;
  29.     _epsilonprec_ = qprecision(q);
  30.     _epsilon_ = qlink(q);
  31.     if (old)
  32.         qfree(old);
  33. }
  34.  
  35.  
  36. /*
  37.  * Return the inverse of one number modulo another.
  38.  * That is, find x such that:
  39.  *    Ax = 1 (mod B)
  40.  * Returns zero if the numbers are not relatively prime (temporary hack).
  41.  */
  42. NUMBER *
  43. qminv(q1, q2)
  44.     NUMBER *q1, *q2;
  45. {
  46.     NUMBER *r;
  47.  
  48.     if (qisfrac(q1) || qisfrac(q2))
  49.         math_error("Non-integers for minv");
  50.     r = qalloc();
  51.     if (zmodinv(q1->num, q2->num, &r->num)) {
  52.         qfree(r);
  53.         return qlink(&_qzero_);
  54.     }
  55.     return r;
  56. }
  57.  
  58.  
  59. /*
  60.  * Return the residue modulo an integer (q3) of an integer (q1) 
  61.  * raised to a positive integer (q2) power.  
  62.  */
  63. NUMBER *
  64. qpowermod(q1, q2, q3)
  65.     NUMBER *q1, *q2, *q3;
  66. {
  67.     NUMBER *r, *t;
  68.     BOOL s;
  69.  
  70.     if (qisfrac(q1) || qisfrac(q2) || qisfrac(q3)) 
  71.         math_error("Non-integers for pmod");
  72.     if (qisneg(q2))
  73.         math_error("Negative power for pmod");
  74.     if (qiszero(q3)) return qpowi(q1, q2);
  75.     s = q3->num.sign;
  76.     q3->num.sign = 0;
  77.     r = qalloc();
  78.     zpowermod(q1->num, q2->num, q3->num, &r->num);
  79.     if (!s || qiszero(r)) return r;
  80.     q3->num.sign = 1;
  81.     t = qadd(r, q3);
  82.     qfree(r);
  83.     return t;
  84. }
  85.  
  86.  
  87. /*
  88.  * Return the power function of one number with another.
  89.  * The power must be integral.
  90.  *    q3 = qpowi(q1, q2);
  91.  */
  92. NUMBER *
  93. qpowi(q1, q2)
  94.     NUMBER *q1, *q2;
  95. {
  96.     register NUMBER *r;
  97.     BOOL invert, sign;
  98.     ZVALUE num, den, z2;
  99.  
  100.     if (qisfrac(q2))
  101.         math_error("Raising number to fractional power");
  102.     num = q1->num;
  103.     den = q1->den;
  104.     z2 = q2->num;
  105.     sign = (num.sign && zisodd(z2));
  106.     invert = z2.sign;
  107.     num.sign = 0;
  108.     z2.sign = 0;
  109.     /*
  110.     * Check for trivial cases first.
  111.     */
  112.     if (ziszero(num) && !ziszero(z2)) {    /* zero raised to a power */
  113.         if (invert)
  114.             math_error("Zero raised to negative power");
  115.         return qlink(&_qzero_);
  116.     }
  117.     if (zisunit(num) && zisunit(den)) {    /* 1 or -1 raised to a power */
  118.         r = (sign ? q1 : &_qone_);
  119.         r->links++;
  120.         return r;
  121.     }
  122.     if (ziszero(z2))    /* raising to zeroth power */
  123.         return qlink(&_qone_);
  124.     if (zisunit(z2)) {    /* raising to power 1 or -1 */
  125.         if (invert)
  126.             return qinv(q1);
  127.         return qlink(q1);
  128.     }
  129.     /*
  130.      * Not a trivial case.  Do the real work.
  131.      */
  132.     r = qalloc();
  133.     if (!zisunit(num))
  134.         zpowi(num, z2, &r->num);
  135.     if (!zisunit(den))
  136.         zpowi(den, z2, &r->den);
  137.     if (invert) {
  138.         z2 = r->num;
  139.         r->num = r->den;
  140.         r->den = z2;
  141.     }
  142.     r->num.sign = sign;
  143.     return r;
  144. }
  145.  
  146.  
  147. /*
  148.  * Given the legs of a right triangle, compute its hypothenuse within
  149.  * the specified error.  This is sqrt(a^2 + b^2).
  150.  */
  151. NUMBER *
  152. qhypot(q1, q2, epsilon)
  153.     NUMBER *q1, *q2, *epsilon;
  154. {
  155.     NUMBER *tmp1, *tmp2, *tmp3;
  156.  
  157.     if (qisneg(epsilon) || qiszero(epsilon))
  158.         math_error("Bad epsilon value for hypot");
  159.     if (qiszero(q1))
  160.         return qabs(q2);
  161.     if (qiszero(q2))
  162.         return qabs(q1);
  163.     tmp1 = qsquare(q1);
  164.     tmp2 = qsquare(q2);
  165.     tmp3 = qadd(tmp1, tmp2);
  166.     qfree(tmp1);
  167.     qfree(tmp2);
  168.     tmp1 = qsqrt(tmp3, epsilon);
  169.     qfree(tmp3);
  170.     return tmp1;
  171. }
  172.  
  173.  
  174. /*
  175.  * Given one leg of a right triangle with unit hypothenuse, calculate
  176.  * the other leg within the specified error.  This is sqrt(1 - a^2).
  177.  * If the wantneg flag is nonzero, then negative square root is returned.
  178.  */
  179. NUMBER *
  180. qlegtoleg(q, epsilon, wantneg)
  181.     NUMBER *q, *epsilon;
  182.     BOOL wantneg;
  183. {
  184.     NUMBER *res, *qtmp1, *qtmp2;
  185.     ZVALUE num;
  186.  
  187.     if (qisneg(epsilon) || qiszero(epsilon))
  188.         math_error("Bad epsilon value for legtoleg");
  189.     if (qisunit(q))
  190.         return qlink(&_qzero_);
  191.     if (qiszero(q)) {
  192.         if (wantneg)
  193.             return qlink(&_qnegone_);
  194.         return qlink(&_qone_);
  195.     }
  196.     num = q->num;
  197.     num.sign = 0;
  198.     if (zrel(num, q->den) >= 0)
  199.         math_error("Leg too large in legtoleg");
  200.     qtmp1 = qsquare(q);
  201.     qtmp2 = qsub(&_qone_, qtmp1);
  202.     qfree(qtmp1);
  203.     res = qsqrt(qtmp2, epsilon);
  204.     qfree(qtmp2);
  205.     if (wantneg) {
  206.         qtmp1 = qneg(res);
  207.         qfree(res);
  208.         res = qtmp1;
  209.     }
  210.     return res;
  211. }
  212.  
  213. /*
  214.  * Compute the square root of a number to within the specified error.
  215.  * If the number is an exact square, the exact result is returned.
  216.  *    q3 = qsqrt(q1, q2);
  217.  */
  218. NUMBER *
  219. qsqrt(q1, epsilon)
  220.     register NUMBER *q1, *epsilon;
  221. {
  222.     long bits, bits2;    /* number of bits of precision */
  223.     int exact;
  224.     NUMBER *r;
  225.     ZVALUE t1, t2;
  226.  
  227.     if (qisneg(q1))
  228.         math_error("Square root of negative number");
  229.     if (qisneg(epsilon) || qiszero(epsilon))
  230.         math_error("Bad epsilon value for sqrt");
  231.     if (qiszero(q1))
  232.         return qlink(&_qzero_);
  233.     if (qisunit(q1))
  234.         return qlink(&_qone_);
  235.     if (qiszero(epsilon) && qisint(q1) && zistiny(q1->num) && (*q1->num.v <= 3))
  236.         return qlink(&_qone_);
  237.     bits = zhighbit(epsilon->den) - zhighbit(epsilon->num) + 1;
  238.     if (bits < 0)
  239.         bits = 0;
  240.     bits2 = zhighbit(q1->num) - zhighbit(q1->den) + 1;
  241.     if (bits2 > 0)
  242.         bits += bits2;
  243.     r = qalloc();
  244.     zshift(q1->num, bits * 2, &t2);
  245.     zmul(q1->den, t2, &t1);
  246.     zfree(t2);
  247.     exact = zsqrt(t1, &t2);
  248.     zfree(t1);
  249.     if (exact) {
  250.         zshift(q1->den, bits, &t1);
  251.         zreduce(t2, t1, &r->num, &r->den);
  252.     } else {
  253.         zquo(t2, q1->den, &t1);
  254.         zfree(t2);
  255.         zbitvalue(bits, &t2);
  256.         zreduce(t1, t2, &r->num, &r->den);
  257.     }
  258.     zfree(t1);
  259.     zfree(t2);
  260.     if (qiszero(r)) {
  261.         qfree(r);
  262.         r = qlink(&_qzero_);
  263.     }
  264.     return r;
  265. }
  266.  
  267.  
  268. /*
  269.  * Calculate the integral part of the square root of a number.
  270.  * Example:  qisqrt(13) = 3.
  271.  */
  272. NUMBER *
  273. qisqrt(q)
  274.     NUMBER *q;
  275. {
  276.     NUMBER *r;
  277.     ZVALUE tmp;
  278.  
  279.     if (qisneg(q))
  280.         math_error("Square root of negative number");
  281.     if (qiszero(q))
  282.         return qlink(&_qzero_);
  283.     if (qisint(q) && zistiny(q->num) && (z1tol(q->num) < 4))
  284.         return qlink(&_qone_);
  285.     r = qalloc();
  286.     if (qisint(q)) {
  287.         (void) zsqrt(q->num, &r->num);
  288.         return r;
  289.     }
  290.     zquo(q->num, q->den, &tmp);
  291.     (void) zsqrt(tmp, &r->num);
  292.     zfree(tmp);
  293.     return r;
  294. }
  295.  
  296.  
  297. /*
  298.  * Return whether or not a number is an exact square.
  299.  */
  300. BOOL
  301. qissquare(q)
  302.     NUMBER *q;
  303. {
  304.     BOOL flag;
  305.  
  306.     flag = zissquare(q->num);
  307.     if (qisint(q) || !flag)
  308.         return flag;
  309.     return zissquare(q->den);
  310. }
  311.  
  312.  
  313. /*
  314.  * Compute the greatest integer of the Kth root of a number.
  315.  * Example:  qiroot(85, 3) = 4.
  316.  */
  317. NUMBER *
  318. qiroot(q1, q2)
  319.     register NUMBER *q1, *q2;
  320. {
  321.     NUMBER *r;
  322.     ZVALUE tmp;
  323.  
  324.     if (qisneg(q2) || qiszero(q2) || qisfrac(q2))
  325.         math_error("Taking number to bad root value");
  326.     if (qiszero(q1))
  327.         return qlink(&_qzero_);
  328.     if (qisone(q1) || qisone(q2))
  329.         return qlink(q1);
  330.     if (qistwo(q2))
  331.         return qisqrt(q1);
  332.     r = qalloc();
  333.     if (qisint(q1)) {
  334.         zroot(q1->num, q2->num, &r->num);
  335.         return r;
  336.     }
  337.     zquo(q1->num, q1->den, &tmp);
  338.     zroot(tmp, q2->num, &r->num);
  339.     zfree(tmp);
  340.     return r;
  341. }
  342.  
  343.  
  344. /*
  345.  * Return the greatest integer of the base 2 log of a number.
  346.  * This is the number such that  1 <= x / log2(x) < 2.
  347.  * Examples:  qilog2(8) = 3, qilog2(1.3) = 1, qilog2(1/7) = -3.
  348.  */
  349. long
  350. qilog2(q)
  351.     NUMBER *q;        /* number to take log of */
  352. {
  353.     long n;            /* power of two */
  354.     int c;            /* result of comparison */
  355.     ZVALUE tmp;        /* temporary value */
  356.  
  357.     if (qisneg(q) || qiszero(q))
  358.         math_error("Non-positive number for log2");
  359.     if (qisint(q))
  360.         return zhighbit(q->num);
  361.     n = zhighbit(q->num) - zhighbit(q->den);
  362.     if (n == 0)
  363.         c = zrel(q->num, q->den);
  364.     else if (n > 0) {
  365.         zshift(q->den, n, &tmp);
  366.         c = zrel(q->num, tmp);
  367.     } else {
  368.         zshift(q->num, -n, &tmp);
  369.         c = zrel(tmp, q->den);
  370.     }
  371.     if (n)
  372.         zfree(tmp);
  373.     if (c < 0)
  374.         n--;
  375.     return n;
  376. }
  377.  
  378.  
  379. /*
  380.  * Return the greatest integer of the base 10 log of a number.
  381.  * This is the number such that  1 <= x / log10(x) < 10.
  382.  * Examples:  qilog10(100) = 2, qilog10(12.3) = 1, qilog10(.023) = -2.
  383.  */
  384. long
  385. qilog10(q)
  386.     NUMBER *q;        /* number to take log of */
  387. {
  388.     long n;            /* log value */
  389.     ZVALUE temp;        /* temporary value */
  390.  
  391.     if (qisneg(q) || qiszero(q))
  392.         math_error("Non-positive number for log10");
  393.     if (qisint(q))
  394.         return zlog10(q->num);
  395.     /*
  396.      * The number is not an integer.
  397.      * Compute the result if the number is greater than one.
  398.      */
  399.     if ((q->num.len > q->den.len) ||
  400.         ((q->num.len == q->den.len) && (zrel(q->num, q->den) > 0))) {
  401.             zquo(q->num, q->den, &temp);
  402.             n = zlog10(temp);
  403.             zfree(temp);
  404.             return n;
  405.     }
  406.     /*
  407.      * Here if the number is less than one.
  408.      * If the number is the inverse of a power of ten, then the obvious answer
  409.      * will be off by one.  Subtracting one if the number is the inverse of an
  410.      * integer will fix it.
  411.      */
  412.     if (zisunit(q->num))
  413.         zsub(q->den, _one_, &temp);
  414.     else
  415.         zquo(q->den, q->num, &temp);
  416.     n = -zlog10(temp) - 1;
  417.     zfree(temp);
  418.     return n;
  419. }
  420.  
  421.  
  422. /*
  423.  * Return the number of digits in a number, ignoring the sign.
  424.  * For fractions, this is the number of digits of its greatest integer.
  425.  * Examples: qdigits(3456) = 4, qdigits(-23.45) = 2, qdigits(.0120) = 1.
  426.  */
  427. long
  428. qdigits(q)
  429.     NUMBER *q;        /* number to count digits of */
  430. {
  431.     long n;            /* number of digits */
  432.     ZVALUE temp;        /* temporary value */
  433.  
  434.     if (qisint(q))
  435.         return zdigits(q->num);
  436.     zquo(q->num, q->den, &temp);
  437.     n = zdigits(temp);
  438.     zfree(temp);
  439.     return n;
  440. }
  441.  
  442.  
  443. /*
  444.  * Return the digit at the specified decimal place of a number represented
  445.  * in floating point.  The lowest digit of the integral part of a number
  446.  * is the zeroth decimal place.  Negative decimal places indicate digits
  447.  * to the right of the decimal point.  Examples: qdigit(1234.5678, 1) = 3,
  448.  * qdigit(1234.5678, -3) = 7.
  449.  */
  450. FLAG
  451. qdigit(q, n)
  452.     NUMBER *q;
  453.     long n;
  454. {
  455.     ZVALUE tenpow, tmp1, tmp2;
  456.     FLAG res;
  457.  
  458.     /*
  459.      * Zero number or negative decimal place of integer is trivial.
  460.      */
  461.     if (qiszero(q) || (qisint(q) && (n < 0)))
  462.         return 0;
  463.     /*
  464.      * For non-negative decimal places, answer is easy.
  465.      */
  466.     if (n >= 0) {
  467.         if (qisint(q))
  468.             return zdigit(q->num, n);
  469.         zquo(q->num, q->den, &tmp1);
  470.         res = zdigit(tmp1, n);
  471.         zfree(tmp1);
  472.         return res;
  473.     }
  474.     /*
  475.      * Fractional value and want negative digit, must work harder.
  476.      */
  477.     ztenpow(-n, &tenpow);
  478.     zmul(q->num, tenpow, &tmp1);
  479.     zfree(tenpow);
  480.     zquo(tmp1, q->den, &tmp2);
  481.     res = zmodi(tmp2, 10L);
  482.     zfree(tmp1);
  483.     zfree(tmp2);
  484.     return res;
  485. }
  486.  
  487.  
  488. /*
  489.  * Return whether or not a bit is set at the specified bit position in a
  490.  * number.  The lowest bit of the integral part of a number is the zeroth
  491.  * bit position.  Negative bit positions indicate bits to the right of the
  492.  * binary decimal point.  Examples: qdigit(17.1, 0) = 1, qdigit(17.1, -1) = 0.
  493.  */
  494. BOOL
  495. qisset(q, n)
  496.     NUMBER *q;
  497.     long n;
  498. {
  499.     NUMBER *qtmp1, *qtmp2;
  500.     ZVALUE ztmp;
  501.     BOOL res;
  502.  
  503.     /*
  504.      * Zero number or negative bit position place of integer is trivial.
  505.      */
  506.     if (qiszero(q) || (qisint(q) && (n < 0)))
  507.         return FALSE;
  508.     /*
  509.      * For non-negative bit positions, answer is easy.
  510.      */
  511.     if (n >= 0) {
  512.         if (qisint(q))
  513.             return zisset(q->num, n);
  514.         zquo(q->num, q->den, &ztmp);
  515.         res = zisset(ztmp, n);
  516.         zfree(ztmp);
  517.         return res;
  518.     }
  519.     /*
  520.      * Fractional value and want negative bit position, must work harder.
  521.      */
  522.     qtmp1 = qscale(q, -n);
  523.     qtmp2 = qint(qtmp1);
  524.     qfree(qtmp1);
  525.     res = ((qtmp2->num.v[0] & 0x01) != 0);
  526.     qfree(qtmp2);
  527.     return res;
  528. }
  529.  
  530.  
  531. /*
  532.  * Compute the factorial of an integer.
  533.  *    q2 = qfact(q1);
  534.  */
  535. NUMBER *
  536. qfact(q)
  537.     register NUMBER *q;
  538. {
  539.     register NUMBER *r;
  540.  
  541.     if (qisfrac(q))
  542.         math_error("Non-integral factorial");
  543.     if (qiszero(q) || zisone(q->num))
  544.         return qlink(&_qone_);
  545.     r = qalloc();
  546.     zfact(q->num, &r->num);
  547.     return r;
  548. }
  549.  
  550.  
  551. /*
  552.  * Compute the product of the primes less than or equal to a number.
  553.  *    q2 = qpfact(q1);
  554.  */
  555. NUMBER *
  556. qpfact(q)
  557.     register NUMBER *q;
  558. {
  559.     NUMBER *r;
  560.  
  561.     if (qisfrac(q))
  562.         math_error("Non-integral factorial");
  563.     r = qalloc();
  564.     zpfact(q->num, &r->num);
  565.     return r;
  566. }
  567.  
  568.  
  569. /*
  570.  * Compute the lcm of all the numbers less than or equal to a number.
  571.  *    q2 = qlcmfact(q1);
  572.  */
  573. NUMBER *
  574. qlcmfact(q)
  575.     register NUMBER *q;
  576. {
  577.     NUMBER *r;
  578.  
  579.     if (qisfrac(q))
  580.         math_error("Non-integral lcmfact");
  581.     r = qalloc();
  582.     zlcmfact(q->num, &r->num);
  583.     return r;
  584. }
  585.  
  586.  
  587. /*
  588.  * Compute the permutation function  M! / (M - N)!.
  589.  */
  590. NUMBER *
  591. qperm(q1, q2)
  592.     register NUMBER *q1, *q2;
  593. {
  594.     register NUMBER *r;
  595.  
  596.     if (qisfrac(q1) || qisfrac(q2))
  597.         math_error("Non-integral arguments for permutation");
  598.     r = qalloc();
  599.     zperm(q1->num, q2->num, &r->num);
  600.     return r;
  601. }
  602.  
  603.  
  604. /*
  605.  * Compute the combinatorial function  M! / (N! * (M - N)!).
  606.  */
  607. NUMBER *
  608. qcomb(q1, q2)
  609.     register NUMBER *q1, *q2;
  610. {
  611.     register NUMBER *r;
  612.  
  613.     if (qisfrac(q1) || qisfrac(q2))
  614.         math_error("Non-integral arguments for combinatorial");
  615.     r = qalloc();
  616.     zcomb(q1->num, q2->num, &r->num);
  617.     return r;
  618. }
  619.  
  620.  
  621. /*
  622.  * Compute the Jacobi function (a / b).
  623.  * -1 => a is not quadratic residue mod b
  624.  * 1 => b is composite, or a is quad residue of b
  625.  * 0 => b is even or b < 0
  626.  */
  627. NUMBER *
  628. qjacobi(q1, q2)
  629.     register NUMBER *q1, *q2;
  630. {
  631.     if (qisfrac(q1) || qisfrac(q2))
  632.         math_error("Non-integral arguments for jacobi");
  633.     return itoq((long) zjacobi(q1->num, q2->num));
  634. }
  635.  
  636.  
  637. /*
  638.  * Compute the Fibonacci number F(n).
  639.  */
  640. NUMBER *
  641. qfib(q)
  642.     register NUMBER *q;
  643. {
  644.     register NUMBER *r;
  645.  
  646.     if (qisfrac(q))
  647.         math_error("Non-integral Fibonacci number");
  648.     r = qalloc();
  649.     zfib(q->num, &r->num);
  650.     return r;
  651. }
  652.  
  653.  
  654. /*
  655.  * Truncate a number to the specified number of decimal places.
  656.  * Specifying zero places makes the result identical to qint.
  657.  * Example: qtrunc(2/3, 3) = .666
  658.  */
  659. NUMBER *
  660. qtrunc(q1, q2)
  661.     NUMBER *q1, *q2;
  662. {
  663.     long places;
  664.     NUMBER *r;
  665.     ZVALUE tenpow, tmp1, tmp2;
  666.  
  667.     if (qisfrac(q2) || qisneg(q2) || !zistiny(q2->num))
  668.         math_error("Bad number of places for qtrunc");
  669.     if (qisint(q1))
  670.         return qlink(q1);
  671.     r = qalloc();
  672.     places = z1tol(q2->num);
  673.     /*
  674.      * Ok, produce the result.
  675.      * First see if we want no places, in which case just take integer part.
  676.      */
  677.     if (places == 0) {
  678.         zquo(q1->num, q1->den, &r->num);
  679.         return r;
  680.     }
  681.     ztenpow(places, &tenpow);
  682.     zmul(q1->num, tenpow, &tmp1);
  683.     zquo(tmp1, q1->den, &tmp2);
  684.     zfree(tmp1);
  685.     if (ziszero(tmp2)) {
  686.         zfree(tmp2);
  687.         return qlink(&_qzero_);
  688.     }
  689.     /*
  690.      * Now reduce the result to the lowest common denominator.
  691.      */
  692.     zgcd(tmp2, tenpow, &tmp1);
  693.     if (zisunit(tmp1)) {
  694.         zfree(tmp1);
  695.         r->num = tmp2;
  696.         r->den = tenpow;
  697.         return r;
  698.     }
  699.     zquo(tmp2, tmp1, &r->num);
  700.     zquo(tenpow, tmp1, &r->den);
  701.     zfree(tmp1);
  702.     zfree(tmp2);
  703.     zfree(tenpow);
  704.     return r;
  705. }
  706.  
  707.  
  708. /*
  709.  * Round a number to the specified number of decimal places.
  710.  * Zero decimal places means round to the nearest integer.
  711.  * Example: qround(2/3, 3) = .667
  712.  */
  713. NUMBER *
  714. qround(q, places)
  715.     NUMBER *q;        /* number to be rounded */
  716.     long places;        /* number of decimal places to round to */
  717. {
  718.     NUMBER *r;
  719.     ZVALUE tenpow, roundval, tmp1, tmp2;
  720.  
  721.     if (places < 0)
  722.         math_error("Negative places for qround");
  723.     if (qisint(q))
  724.         return qlink(q);
  725.     /*
  726.      * Calculate one half of the denominator, ignoring fractional results.
  727.      * This is the value we will add in order to cause rounding.
  728.      */
  729.     zshift(q->den, -1L, &roundval);
  730.     roundval.sign = q->num.sign;
  731.     /*
  732.      * Ok, now do the actual work to produce the result.
  733.      */
  734.     r = qalloc();
  735.     ztenpow(places, &tenpow);
  736.     zmul(q->num, tenpow, &tmp2);
  737.     zadd(tmp2, roundval, &tmp1);
  738.     zfree(tmp2);
  739.     zfree(roundval);
  740.     zquo(tmp1, q->den, &tmp2);
  741.     zfree(tmp1);
  742.     if (ziszero(tmp2)) {
  743.         zfree(tmp2);
  744.         return qlink(&_qzero_);
  745.     }
  746.     /*
  747.      * Now reduce the result to the lowest common denominator.
  748.      */
  749.     zgcd(tmp2, tenpow, &tmp1);
  750.     if (zisunit(tmp1)) {
  751.         zfree(tmp1);
  752.         r->num = tmp2;
  753.         r->den = tenpow;
  754.         return r;
  755.     }
  756.     zquo(tmp2, tmp1, &r->num);
  757.     zquo(tenpow, tmp1, &r->den);
  758.     zfree(tmp1);
  759.     zfree(tmp2);
  760.     zfree(tenpow);
  761.     return r;
  762. }
  763.  
  764.  
  765. /*
  766.  * Truncate a number to the specified number of binary places.
  767.  * Specifying zero places makes the result identical to qint.
  768.  */
  769. NUMBER *
  770. qbtrunc(q1, q2)
  771.     NUMBER *q1, *q2;
  772. {
  773.     long places, twopow;
  774.     NUMBER *r;
  775.     ZVALUE tmp1, tmp2;
  776.  
  777.     if (qisfrac(q2) || qisneg(q2) || !zistiny(q2->num))
  778.         math_error("Bad number of places for qtrunc");
  779.     if (qisint(q1))
  780.         return qlink(q1);
  781.     r = qalloc();
  782.     places = z1tol(q2->num);
  783.     /*
  784.      * Ok, produce the result.
  785.      * First see if we want no places, in which case just take integer part.
  786.      */
  787.     if (places == 0) {
  788.         zquo(q1->num, q1->den, &r->num);
  789.         return r;
  790.     }
  791.     zshift(q1->num, places, &tmp1);
  792.     zquo(tmp1, q1->den, &tmp2);
  793.     zfree(tmp1);
  794.     if (ziszero(tmp2)) {
  795.         zfree(tmp2);
  796.         return qlink(&_qzero_);
  797.     }
  798.     /*
  799.      * Now reduce the result to the lowest common denominator.
  800.      */
  801.     if (zisodd(tmp2)) {
  802.         r->num = tmp2;
  803.         zbitvalue(places, &r->den);
  804.         return r;
  805.     }
  806.     twopow = zlowbit(tmp2);
  807.     if (twopow > places)
  808.         twopow = places;
  809.     places -= twopow;
  810.     zshift(tmp2, -twopow, &r->num);
  811.     zfree(tmp2);
  812.     zbitvalue(places, &r->den);
  813.     return r;
  814. }
  815.  
  816.  
  817. /*
  818.  * Round a number to the specified number of binary places.
  819.  * Zero binary places means round to the nearest integer.
  820.  */
  821. NUMBER *
  822. qbround(q, places)
  823.     NUMBER *q;        /* number to be rounded */
  824.     long places;        /* number of binary places to round to */
  825. {
  826.     long twopow;
  827.     NUMBER *r;
  828.     ZVALUE roundval, tmp1, tmp2;
  829.  
  830.     if (places < 0)
  831.         math_error("Negative places for qbround");
  832.     if (qisint(q))
  833.         return qlink(q);
  834.     r = qalloc();
  835.     /*
  836.      * Calculate one half of the denominator, ignoring fractional results.
  837.      * This is the value we will add in order to cause rounding.
  838.      */
  839.     zshift(q->den, -1L, &roundval);
  840.     roundval.sign = q->num.sign;
  841.     /*
  842.      * Ok, produce the result.
  843.      */
  844.     zshift(q->num, places, &tmp1);
  845.     zadd(tmp1, roundval, &tmp2);
  846.     zfree(roundval);
  847.     zfree(tmp1);
  848.     zquo(tmp2, q->den, &tmp1);
  849.     zfree(tmp2);
  850.     if (ziszero(tmp1)) {
  851.         zfree(tmp1);
  852.         return qlink(&_qzero_);
  853.     }
  854.     /*
  855.      * Now reduce the result to the lowest common denominator.
  856.      */
  857.     if (zisodd(tmp1)) {
  858.         r->num = tmp1;
  859.         zbitvalue(places, &r->den);
  860.         return r;
  861.     }
  862.     twopow = zlowbit(tmp1);
  863.     if (twopow > places)
  864.         twopow = places;
  865.     places -= twopow;
  866.     zshift(tmp1, -twopow, &r->num);
  867.     zfree(tmp1);
  868.     zbitvalue(places, &r->den);
  869.     return r;
  870. }
  871.  
  872.  
  873. /*
  874.  * Approximate a number by using binary rounding with the minimum number
  875.  * of binary places so that the resulting number differs from the
  876.  * original number by no more than epsilon/2.
  877.  */
  878.  
  879. NUMBER *
  880. qbappr(q, e)
  881.     NUMBER *q, *e;
  882. {
  883.     long bits;
  884.  
  885.     if (qiszero(e))
  886.         return qlink(q);
  887.     if (qisneg(e))
  888.         math_error("Bad epsilon value for qbappr");
  889.     if (e == _epsilon_)
  890.         bits = _epsilonprec_;
  891.     else
  892.         bits = qprecision(e);
  893.     return qbround(q, bits);
  894. }
  895.  
  896.  
  897. /*
  898.  * Approximate a number using continued fractions until the approximation
  899.  * error is less than the specified value.  If a NULL pointer is given
  900.  * for the error value, then the closest simpler fraction is returned.
  901.  */
  902. NUMBER *
  903. qcfappr(q, e)
  904.     NUMBER *q, *e;
  905. {
  906.     NUMBER qtest, *qtmp;
  907.     ZVALUE u1, u2, u3, v1, v2, v3, t1, t2, t3, qq, tt;
  908.     int i;
  909.     BOOL haveeps;
  910.  
  911.     haveeps = TRUE;
  912.     if (e == NULL) {
  913.         haveeps = FALSE;
  914.         e = &_qzero_;
  915.     }
  916.     if (qisneg(e))
  917.         math_error("Negative epsilon for cfappr");
  918.     if (qisint(q) || zisunit(q->num) || (haveeps && qiszero(e)))
  919.         return qlink(q);
  920.     u1 = _one_;
  921.     u2 = _zero_;
  922.     u3 = q->num;
  923.     u3.sign = 0;
  924.     v1 = _zero_;
  925.     v2 = _one_;
  926.     v3 = q->den;
  927.     while (!ziszero(v3)) {
  928.         if (!ziszero(u2) && !ziszero(u1)) {
  929.             qtest.num = u2;
  930.             qtest.den = u1;
  931.             qtest.den.sign = 0;
  932.             qtest.num.sign = q->num.sign;
  933.             qtmp = qsub(q, &qtest);
  934.             qtest = *qtmp;
  935.             qtest.num.sign = 0;
  936.             i = qrel(&qtest, e);
  937.             qfree(qtmp);
  938.             if (i <= 0)
  939.                 break;
  940.         }
  941.         zquo(u3, v3, &qq);
  942.         zmul(qq, v1, &tt); zsub(u1, tt, &t1); zfree(tt);
  943.         zmul(qq, v2, &tt); zsub(u2, tt, &t2); zfree(tt);
  944.         zmul(qq, v3, &tt); zsub(u3, tt, &t3); zfree(tt);
  945.         zfree(qq); zfree(u1); zfree(u2);
  946.         if ((u3.v != q->num.v) && (u3.v != q->den.v))
  947.             zfree(u3);
  948.         u1 = v1; u2 = v2; u3 = v3;
  949.         v1 = t1; v2 = t2; v3 = t3;
  950.     }
  951.     if (u3.v != q->den.v)
  952.         zfree(u3);
  953.     zfree(v1);
  954.     zfree(v2);
  955.     i = ziszero(v3);
  956.     zfree(v3);
  957.     if (i && haveeps) {
  958.         zfree(u1);
  959.         zfree(u2);
  960.         return qlink(q);
  961.     }
  962.     qtest.num = u2;
  963.     qtest.den = u1;
  964.     qtest.den.sign = 0;
  965.     qtest.num.sign = q->num.sign;
  966.     qtmp = qcopy(&qtest);
  967.     zfree(u1);
  968.     zfree(u2);
  969.     return qtmp;
  970. }
  971.  
  972.  
  973. /*
  974.  * Return an indication on whether or not two fractions are approximately
  975.  * equal within the specified epsilon. Returns negative if the absolute value
  976.  * of the difference between the two numbers is less than epsilon, zero if
  977.  * the difference is equal to epsilon, and positive if the difference is
  978.  * greater than epsilon.
  979.  */
  980. FLAG
  981. qnear(q1, q2, epsilon)
  982.     NUMBER *q1, *q2, *epsilon;
  983. {
  984.     int res;
  985.     NUMBER qtemp, *qq;
  986.  
  987.     if (qisneg(epsilon))
  988.         math_error("Negative epsilon for near");
  989.     if (q1 == q2) {
  990.         if (qiszero(epsilon))
  991.             return 0;
  992.         return -1;
  993.     }
  994.     if (qiszero(epsilon))
  995.         return qcmp(q1, q2);
  996.     if (qiszero(q2)) {
  997.         qtemp = *q1;
  998.         qtemp.num.sign = 0;
  999.         return qrel(&qtemp, epsilon);
  1000.     }
  1001.     if (qiszero(q1)) {
  1002.         qtemp = *q2;
  1003.         qtemp.num.sign = 0;
  1004.         return qrel(&qtemp, epsilon);
  1005.     }
  1006.     qq = qsub(q1, q2);
  1007.     qtemp = *qq;
  1008.     qtemp.num.sign = 0;
  1009.     res = qrel(&qtemp, epsilon);
  1010.     qfree(qq);
  1011.     return res;
  1012. }
  1013.  
  1014.  
  1015. /*
  1016.  * Compute the gcd (greatest common divisor) of two numbers.
  1017.  *    q3 = qgcd(q1, q2);
  1018.  */
  1019. NUMBER *
  1020. qgcd(q1, q2)
  1021.     register NUMBER *q1, *q2;
  1022. {
  1023.     ZVALUE z;
  1024.     NUMBER *q;
  1025.  
  1026.     if (q1 == q2)
  1027.         return qabs(q1);
  1028.     if (qisfrac(q1) || qisfrac(q2)) {
  1029.         q = qalloc();
  1030.         zgcd(q1->num, q2->num, &q->num);
  1031.         zlcm(q1->den, q2->den, &q->den);
  1032.         return q;
  1033.     }
  1034.     if (qiszero(q1))
  1035.         return qabs(q2);
  1036.     if (qiszero(q2))
  1037.         return qabs(q1);
  1038.     if (qisunit(q1) || qisunit(q2))
  1039.         return qlink(&_qone_);
  1040.     zgcd(q1->num, q2->num, &z);
  1041.     if (zisunit(z)) {
  1042.         zfree(z);
  1043.         return qlink(&_qone_);
  1044.     }
  1045.     q = qalloc();
  1046.     q->num = z;
  1047.     return q;
  1048. }
  1049.  
  1050.  
  1051. /*
  1052.  * Compute the lcm (least common multiple) of two numbers.
  1053.  *    q3 = qlcm(q1, q2);
  1054.  */
  1055. NUMBER *
  1056. qlcm(q1, q2)
  1057.     register NUMBER *q1, *q2;
  1058. {
  1059.     NUMBER *q;
  1060.  
  1061.     if (qiszero(q1) || qiszero(q2))
  1062.         return qlink(&_qzero_);
  1063.     if (q1 == q2)
  1064.         return qabs(q1);
  1065.     if (qisunit(q1))
  1066.         return qabs(q2);
  1067.     if (qisunit(q2))
  1068.         return qabs(q1);
  1069.     q = qalloc();
  1070.     zlcm(q1->num, q2->num, &q->num);
  1071.     if (qisfrac(q1) || qisfrac(q2))
  1072.         zgcd(q1->den, q2->den, &q->den);
  1073.     return q;
  1074. }
  1075.  
  1076.  
  1077. /*
  1078.  * Remove all occurances of the specified factor from a number.
  1079.  * Returned number is always positive.
  1080.  */
  1081. NUMBER *
  1082. qfacrem(q1, q2)
  1083.     NUMBER *q1, *q2;
  1084. {
  1085.     long count;
  1086.     ZVALUE tmp;
  1087.     NUMBER *r;
  1088.  
  1089.     if (qisfrac(q1) || qisfrac(q2))
  1090.         math_error("Non-integers for factor removal");
  1091.     count = zfacrem(q1->num, q2->num, &tmp);
  1092.     if (zisunit(tmp)) {
  1093.         zfree(tmp);
  1094.         return qlink(&_qone_);
  1095.     }
  1096.     if (count == 0) {
  1097.         zfree(tmp);
  1098.         return qlink(q1);
  1099.     }
  1100.     r = qalloc();
  1101.     r->num = tmp;
  1102.     return r;
  1103. }
  1104.  
  1105.  
  1106. /*
  1107.  * Divide one number by the gcd of it with another number repeatedly until
  1108.  * the number is relatively prime.
  1109.  */
  1110. NUMBER *
  1111. qgcdrem(q1, q2)
  1112.     NUMBER *q1, *q2;
  1113. {
  1114.     ZVALUE tmp;
  1115.     NUMBER *r;
  1116.  
  1117.     if (qisfrac(q1) || qisfrac(q2))
  1118.         math_error("Non-integers for gcdrem");
  1119.     zgcdrem(q1->num, q2->num, &tmp);
  1120.     if (zisunit(tmp)) {
  1121.         zfree(tmp);
  1122.         return qlink(&_qone_);
  1123.     }
  1124.     if (zcmp(q1->num, tmp) == 0) {
  1125.         zfree(tmp);
  1126.         return qlink(q1);
  1127.     }
  1128.     r = qalloc();
  1129.     r->num = tmp;
  1130.     return r;
  1131. }
  1132.  
  1133.  
  1134. /*
  1135.  * Return the lowest prime factor of a number.
  1136.  * Search is conducted for the specified number of primes.
  1137.  * Returns one if no factor was found.
  1138.  */
  1139. NUMBER *
  1140. qlowfactor(q1, q2)
  1141.     NUMBER *q1, *q2;
  1142. {
  1143.     if (qisfrac(q1) || qisfrac(q2))
  1144.         math_error("Non-integers for lowfactor");
  1145.     return itoq(zlowfactor(q1->num, ztoi(q2->num)));
  1146. }
  1147.  
  1148.  
  1149. /*
  1150.  * Return the number of places after the decimal point needed to exactly
  1151.  * represent the specified number as a real number.  Integers return zero,
  1152.  * and non-terminating decimals return minus one.  Examples:
  1153.  *    qplaces(1/7)=-1, qplaces(3/10)= 1, qplaces(1/8)=3, qplaces(4)=0.
  1154.  */
  1155. long
  1156. qplaces(q)
  1157.     NUMBER *q;
  1158. {
  1159.     long twopow, fivepow;
  1160.     HALF fiveval[2];
  1161.     ZVALUE five;
  1162.     ZVALUE tmp;
  1163.  
  1164.     if (qisint(q))    /* no decimal places if number is integer */
  1165.         return 0;
  1166.     /*
  1167.      * The number of decimal places of a fraction in lowest terms is finite
  1168.      * if an only if the denominator is of the form 2^A * 5^B, and then the
  1169.      * number of decimal places is equal to MAX(A, B).
  1170.      */
  1171.     five.sign = 0;
  1172.     five.len = 1;
  1173.     five.v = fiveval;
  1174.     fiveval[0] = 5;
  1175.     fivepow = zfacrem(q->den, five, &tmp);
  1176.     if (!zisonebit(tmp)) {
  1177.         zfree(tmp);
  1178.         return -1;
  1179.     }
  1180.     twopow = zlowbit(tmp);
  1181.     zfree(tmp);
  1182.     if (twopow < fivepow)
  1183.         twopow = fivepow;
  1184.     return twopow;
  1185. }
  1186.  
  1187.  
  1188. /*
  1189.  * Perform a probabilistic primality test (algorithm P in Knuth).
  1190.  * Returns FALSE if definitely not prime, or TRUE if probably prime.
  1191.  * Second arg determines how many times to check for primality.
  1192.  */
  1193. BOOL
  1194. qprimetest(q1, q2)
  1195.     NUMBER *q1, *q2;
  1196. {
  1197.     if (qisfrac(q1) || qisfrac(q2) || qisneg(q2))
  1198.         math_error("Bad arguments for qprimetest");
  1199.     return zprimetest(q1->num, qtoi(q2));
  1200. }
  1201.  
  1202.  
  1203. /*
  1204.  * Return a trivial hash value for a number.
  1205.  */
  1206. HASH
  1207. qhash(q)
  1208.     NUMBER *q;
  1209. {
  1210.     HASH hash;
  1211.  
  1212.     hash = zhash(q->num);
  1213.     if (qisfrac(q))
  1214.         hash += zhash(q->den) * 2000003;
  1215.     return hash;
  1216. }
  1217.  
  1218. /* END CODE */
  1219.